Journal list menu

ARTICLE
Open Access

Alpine butterflies want to fly high: Species and communities shift upwards faster than their host plants

First published: 13 August 2022
Citations: 1
Handling Editor: Mariano A. Rodriguez-Cabal

Funding information: Bavarian Ministry for Science and the Arts

Abstract

Despite sometimes strong codependencies of insect herbivores and plants, the responses of individual taxa to accelerating climate change are typically studied in isolation. For this reason, biotic interactions that potentially limit species in tracking their preferred climatic niches are ignored. Here, we chose butterflies as a prominent representative of herbivorous insects to investigate the impacts of temperature changes and their larval host plant distributions along a 1.4-km elevational gradient in the German Alps. Following a sampling protocol of 2009, we revisited 33 grassland plots in 2019 over an entire growing season. We quantified changes in butterfly abundance and richness by repeated transect walks on each plot and disentangled the direct and indirect effects of locally assessed temperature, site management, and larval and adult food resource availability on these patterns. Additionally, we determined elevational range shifts of butterflies and host plants at both the community and species level. Comparing the two sampled years (2009 and 2019), we found a severe decline in butterfly abundance and a clear upward shift of butterflies along the elevational gradient. We detected shifts in the peak of species richness, community composition, and at the species level, whereby mountainous species shifted particularly strongly. In contrast, host plants showed barely any change, neither in connection with species richness nor individual species shifts. Further, temperature and host plant richness were the main drivers of butterfly richness, with change in temperature best explaining the change in richness over time. We concluded that host plants were not yet hindering butterfly species and communities from shifting upwards. However, the mismatch between butterfly and host plant shifts might become a problem for this very close plant–herbivore relationship, especially toward higher elevations, if butterflies fail to adapt to new host plants. Further, our results support the value of conserving traditional extensive pasture use as a promoter of host plant and, hence, butterfly richness.

INTRODUCTION

Climate change is globally one of the main drivers restructuring species communities (Halsch et al., 2021; Hill et al., 2021; Parmesan, 2006). Especially in mountainous regions, species communities often show particularly pronounced and rapid reactions to climate change (Bertrand et al., 2011; Lawler et al., 2009; Pompe et al., 2010): On the one hand because mountains harbor an extraordinary diversity of endemic and cold-adapted species (Laiolo et al., 2018; Viterbi et al., 2013), while facing a faster-than-average temperature increase (Pepin et al., 2015), and on the other hand because the steep natural temperature gradient along the mountain slopes offers species communities the opportunity to adapt comparatively quickly to climatic changes by, for example, shifting their ranges to higher elevations (Chen et al., 2009; Lenoir et al., 2008; Molina-Martínez et al., 2016). However, species ranges are not only shaped by temperature (Normand et al., 2009); biotic interactions play a major role in the manifestation of biotic responses to climatic change (Alexander et al., 2015; Cahill et al., 2013). Nevertheless, studies investigating the responses of interacting species on the community level in the context of climate change remain rare.

Ectotherms like insects are sensitive to temperature changes and together with short life cycles and fast generation turnover, insects react faster than other organism groups to changes in temperature (Lenoir et al., 2020; Robinet & Roques, 2010). Therefore, they are a good indicator to track climate change impacts. An excellent representative of insects for monitoring are butterflies, because as terrestrial flying insects, they are quite mobile and were also observed to react faster than other organism groups in the Alps (Vitasse et al., 2021). Further, butterflies have been of public interest for a very long time, so that—in comparison to many other insect groups—many data sets around the world are available, also over long time periods (Settele et al., 2008). However, with their herbivorous lifestyle, they take a special position among insects, insofar as many species are dependent on certain larval host plants. Several studies in different mountain ranges, including the Alps, revealed either the distributional shifts of plant (Lenoir et al., 2008; Roth et al., 2014; Vitasse et al., 2021) or butterfly communities under climate change (Cerrato et al., 2019; Forister et al., 2010; Rödder et al., 2021). Although modeling approaches suggest that host plant availability might limit butterfly distributions and climate change can disrupt trophically interacting species (Hanspach et al., 2014; Schweiger et al., 2008), only few empirical studies consider temperature and host plant diversity as potentially interacting drivers of butterfly diversity (Gutiérrez et al., 2016; Menéndez et al., 2007), and these have rarely examined the temporal dynamics of butterfly and host plant diversity in response to climate change.

Highly protected areas, such as national parks, are ideal for studying the natural responses of host plants and butterflies to climate change. Elsewhere, such responses are typically modified or entirely masked by the negative effects of intensive land use or other anthropogenic activities on natural communities (Guo et al., 2018; Parks et al., 2020). Further, by focusing on grasslands, the impact by differing habitat types like, for example, different habitat structures and microclimates can be diminished. Grasslands are of high importance as hotspots of biodiversity (Habel et al., 2013). They constitute the main habitat for butterflies in general (WallisDeVries & van Swaay, 2009) and in the Alps (Bräu et al., 2013). However, mountainous grasslands below the timberline are mainly a heritage of extensive anthropogenic land-use practices in Central Europe (Ellenberg & Leuschner, 2010), becoming more and more abandoned due to declining economic viability (Laiolo et al., 2004). Because missing extensive management negatively affects butterfly species richness (Jerrentrup et al., 2014), we here aim to clarify the magnitude of this impact.

In this study, we revisited 33 grassland sites in the National Park Berchtesgaden in the German Alps, to assess butterfly and the respective host plant communities along elevational gradients, one decade after the first assessment (Hoiss et al., 2012; Leingärtner et al., 2014). Our community comparison took place within the currently hottest decade on record (2009 vs. 2019) over which temperature rose around ~0.3°C globally (Li et al., 2021) and even 0.1°C more in the Alps (Nigrelli & Chiarle, 2021). We investigated the following questions: (1) How did butterfly and host plant richness change along the elevational gradient in this extraordinary warm decade? (2) What is the relative importance of temperature, management, and host plant richness as potential drivers of butterfly richness changes in space and time? (3) Do butterfly communities and species shift along the elevational gradient in magnitudes that mirror the changes in temperature, or are they limited by their association with host plants?

MATERIAL AND METHODS

Study area

This study was carried out in the National Park Berchtesgaden and its surroundings, located in the southeasternmost corner of Germany as part of the Berchtesgaden Alps, which belong to the Northern Limestone Alps. The national park is characterized by steep valleys and high mountains ranging from 600 up to 2700 m above sea level (MASL). The mean annual temperature (MAT) averaged over the last 10 years (2011–2020) varies between −0.8 and 9.7°C, depending on the elevation, with an annual precipitation of 1240–2881 mm (extracted from 1 × 1-km grid data of yearly means from the German Meteorological Service [DWD]). The higher elevations above the timberline host mainly alpine meadows, whereas coniferous forests interspersed with pastures dominate lower elevations. Alpine pastoral use has a centuries-long tradition in this region. However, some meadows are abandoned, whereas others are still traditionally managed.

In 2009, 34 grassland sites were selected, covering an elevational range from 600 up to approximately 2000 MASL (Leingärtner et al., 2014), which was the highest elevation at which a closed grass layer can be found. The sites were arranged along five transects with approximately 250 m in altitude in between. Twenty-nine sites were located inside the national park and five in the close surroundings to extend the elevational gradient towards the lowlands (Leingärtner et al., 2014). All grasslands were unused or extensively managed by extensive grazing or one late summer cut without any fertilization. In 2019, 33 of these 34 grassland sites could be revisited (Figure 1c). Owing to the national park status of most sites, management remained stable within the considered decade.

Details are in the caption following the image
(a) Location of study area in the southeasternmost corner of Germany in National Park Berchtesgaden (NP BGD, green) and its surroundings. (b) Temperature differences in degrees Celsius between 2009 and 2019 along the elevational gradient in NP BGD. There was a considerable difference between mean summer temperature (MST) (purple) and mean annual temperature (MAT) (blue) increase. The darker shades show the absolute difference between the two sampling years (2009/2019), and the lighter ones show the gradual temperature change over the considered decade (2009–2019). Temperatures were predicted with generalized additive models (GAMs) from adjacent climate stations for each study site. (c) Climate stations used for temperature modeling were spread throughout the national park and its surroundings (orange). Study sites (light blue) were arranged along five transects (dark blue) plus four valley plots at lower elevations. Gray shades indicate elevation in 200-m intervals, ranging from 400 (black) up to 2800 m above sea level (MASL) (white) (based on a digital elevation model with 1 m grid cell width, provided by national park administration).

Data collection

Butterfly monitoring

We monitored diurnal butterflies (Lepidoptera: Hesperiidae, Lycaenidae, Nymphalidae, Papilionidae, Pieridae, Riodinidae) from 2 May to 16 September 2019, exactly following the sampling protocol from 2009 (8 May to 10 September) (Leingärtner et al., 2014). Sites below 1200 MASL were sampled six times and those above five times due to the shorter vegetation period. This had no effect on the average species detection rate (plots five times sampled 73%–100%, six times sampled 76%–99%, based on detected species richness/abundance-based coverage estimator [ACE]; see section Statistical analyses for further details), indicating a similar coverage of the species community at all sites. On each site (60 × 60 m), we conducted randomized transect walks of 400 m in 32 min, subdivided into 8 × 4 min to calculate species saturation. The considered transect corridor spanned 2.5 m to the right and left of the observer as well as 5 m to the front. We sampled between 9:00 AM and 6:00 PM either when the sun was shining or when temperature reached at least 17°C on cloudy days. Further, sampling was restricted to winds of force ≤4 on the Beaufort scale. Additionally, we estimated flower cover as the share of the whole plot surface covered by flowers for every visit. All sites were visited at regular intervals over the season and in random order. Individuals were caught with a butterfly net, identified on site, and immediately released afterwards. During identification, the timer was stopped. We collected only one voucher specimen per species and not clearly identifiable specimens for further identification. In particular, the genus Erebia needed further examination by a trained taxonomist. Closely related species that we could not clearly differentiate by visual inspection were summarized in species complexes. Those were Aricia agestis and A. artaxerxes, Colias alfacariensis and C. hyale, Leptidea sinapis, L. juvernica and L. reali, and Pyrgus alveus with P. cacaliae and P. serratulae. Further, we did not differentiate between Pieris bryoniae and P. napi. Although females look clearly different, males cannot be clearly distinguished owing to wing bleaching and hybridization. Species identification and nomenclature were based on Settele et al. (2015) and complemented by Paolucci (2013) for alpine species. In total, we recorded 4595 individuals of 65 butterfly species in 2009 and 2449 individuals of 68 species in 2019, whereby 62 species occurred in both years.

Vegetation assessment

Plant species and cover were assessed once between June and July 2019 on all 33 sites of the butterfly monitoring by the same botanists as in 2009 (Hoiss et al., 2012), exactly repeating the sampling protocol. However, one site had to be excluded in the end owing to feeding damage by sheep shortly before data collection in 2019, preventing species identification. We randomly distributed ten 2 × 2-m subplots within each study site at a minimum distance of 5 m between subplots. In each subplot, we assessed all vascular plant species and estimated their respective cover, using the Domin scale, which combines abundance and cover information at a more detailed level than the widely used Braun-Blanquet scale (Mueller-Dombois & Ellenberg, 1974).

For further analyses, we filtered the data for records of larval host plants of butterflies. The plants were defined as such if mentioned in Paolucci (2013) or Reinhardt et al. (2020). In case they only mentioned the name of a genus or family, we considered all associated species as host plants. Associations of butterflies with host plants are presented in Appendix S1: Table S1.

Abiotic factors

For the quantification of long-term temperature changes and between-year temperature differences, we predicted daily mean temperatures between 1 January 2009 and 31 December 2019 for all study sites. We used temperature records from 11 (2009) to 16 (2019) adjacent climate stations, located in the National Park Berchtesgaden and its surroundings, covering an elevational range from 653 to 2645 MASL (Figure 1). For the prediction, temperature, which was recorded in 10-min intervals, was first averaged per day and climate station. We then fitted generalized additive models (GAMs) with a smooth term of elevation, day length in minutes (as a proxy for season), and date as a factor. As the number of climate stations increased over the decade, we modeled temperature data for each year separately and with the maximal number of available stations. Deviance explained by the GAMs ranged between 96 and 99% depending on the year (see also Appendix S1: Figure S1, which demonstrates measured and predicted temperatures in 2019 as an example). Since temperature models are based solely on elevation, day length, and date but do not correct for topographic site properties like shading, slope, or exposure, there might be some site-specific variation that is not depicted in our modeled temperature data. The exact elevation of sites and climate stations was extracted from a digital elevation model with 1-m resolution, provided by the national park administration.

Statistical analyses

We performed all statistical analyses using R version 4.1.1 (R Core Team, 2021). Daily temperatures were averaged per year (mean annual temperature [MAT]) from 2009 to 2019 and for the respective summer months (May–September) (mean summer temperature [MST]). We chose MST as a predictor variable because it had a better fit in our models than MAT, reflecting the sampling period and the main flight time of adult butterflies. In addition to the absolute temperature difference between the two sampling years (2009/2019), we estimated the gradual temperature change over the decade (2009–2019) by extracting the slope of linear models (MAT/MST ~ year) (Appendix S1: Figure S2), showing a steeper trend for the temperature increase over the decade than the absolute difference between the sampling years. The changes in temperature differences along the elevational gradient were analyzed using GAMs with the package mgcv version 1.8–38 (Wood, 2017), fitted with a smoothing function of elevation and setting the basis dimension k of the smoothing parameter to four (Figure 1b, Table 1). By setting k = 4, we limit the number of turns the curve is allowed to make when fitted to the data, equaling a polynomial function of maximum fourth order to avoid overfitting and keep the number of turns in an explainable dimension. We used k = 4 for all GAMs and generalized additive mixed models (GAMMs) described in the following paragraphs. We determined the local temperature lapse rate along the elevational gradient by extracting the slope of a linear model (temperature and elevation), based on mean daily temperatures from nearby climate stations, averaged over the considered decade per station (Appendix S1: Figure S3). The lapse rate and the change of temperature between or across years were used to calculate expected species range shifts and to compare them with observed shifts.

TABLE 1. Summary statistics.
Response variable Predictor EDF Ref.df p-value Family Link function Adjustted R2 Deviance explained
∆MAT 2009/2019 Elevation 2.994 2.994 <0.001*** Gaussian Identity 0.999 99.9%
∆MST 2009/2019 Elevation 2.994 2.994 <0.001*** Gaussian Identity 0.999 99.9%
∆MAT 2009–2019 Elevation 2.997 2.997 <0.001*** Gaussian Identity 0.999 99.9%
∆MST 2009–2019 Elevation 2.997 2.997 <0.001*** Gaussian Identity 0.999 99.9%
Butterfly abundance Elevation 1.111 1.115 0.983 n.s. Poisson Log 0.721 84.4%
Elevation:year 1.000 1.000 <0.001***
Study site (random term) 29.878 31.000 <0.001***
Butterfly species no. Elevation 1.000 1.000 <0.001*** Poisson Log 0.612 72.5%
Elevation:year 2.159 2.475 0.018*
Study site (random term) 16.113 31.000 <0.001***
Host plant species no. Elevation 1.783 1.846 0.007** Poisson Log 0.890 93.1%
Elevation:year 1.000 1.000 0.418 n.s.
Study site (random term) 24.637 30.000 <0.001***
Butterfly community score Elevation 2.460 2.539 <0.001*** Gaussian Identity 0.963 97.9%
Elevation:year 1.865 2.234 <0.001***
Study site (random term) 22.919 31.000 <0.001***
  • Note: These statistics relate to generalized additive models (GAMs) for absolute (2009/2019) and gradual (2009–2019) changes in mean annual temperature (MAT) and mean summer temperature (MST), fitted with a smoothing function of elevation (N = 32) (Figure 1b), as well as generalized additive mixed models (GAMMs) for butterfly abundance, butterfly, and host plant species richness and the community elevation score, each fitted with a smoothing function of elevation and the interaction of elevation and year (2009, 2019), including study site as a random term (N = 64). Significance codes: 0 ≤ “***” < 0.001 < “**” < 0.01 < “*” < 0.05 < “.” < 0.1 < “n.s.” < 1.

We pooled butterfly species records per site for each year and averaged the estimated flower cover per site over all visits of the year. The scale values of the vegetation survey for each host plant species were transformed into an averaged percentage cover value of all 10 subplots per site, using the Domin 2.6 transformation from Currall (1987). Species richness of butterflies and host plants was calculated with the vegan package version 2.5–7 (Oksanen et al., 2020) and analyzed using GAMMs with the mgcv package version 1.8–38 (Wood, 2017). The models were fitted with a smoothing function of elevation, year as a factor (2009, 2019), and the interaction of elevation and year, including study site as a random intercept effect and setting data family to Poisson. We used the observed species richness of butterflies for calculations because it was highly correlated (r = 0.98) with the abundance-based coverage estimator (Chao & Lee, 1992) (ACE, calculated with the package fossil V. 0.4.0 [Vavrek, 2011]) showing a high species richness saturation (73%–100%). Thus, the data seemed to reflect well the local species community. Because of the many “rare” species with few individuals, especially at high elevations with harsh weather conditions, ACE was most appropriate for validation in comparison to other indices like Chao1. It recognizes those “rare” species as such (Chao & Shen, 2004) instead of overestimating species richness at plots with many low numbers of individuals.

We used piecewise structural equation models (SEMs), provided by the piecewiseSEM package version 2.1.2 (Lefcheck, 2016), to identify (a) drivers of butterfly species richness and (b) drivers of the change in richness between years. In comparison to traditional SEMs, these models allow a wider variety of model structures and can account for various data distributions (e.g., Poisson) or relatively small sample sizes as typical issues in ecological research. We also considered the latter in model selection using the Akaike information criterion with a correction for small sample sizes (AICc) for decision-making. We selected four potential drivers of butterfly richness: Mean summer temperature (MST), management, adult food resources (flower cover), and number of larval food plants (host plant richness). We chose MST because the model had a higher explanatory value and better fit (lower AICc) than with MAT. Management per site was classified either as extensively used (1), including one late summer cut and extensive grazing, or unused (0). In addition to host plant richness, host plant cover was included initially but disregarded in the end because it added no explanatory value to the model. For the second SEM, explaining the change in butterfly richness, we included the absolute differences of all explaining variables in the SEM (e.g., ∆MST2019–2009). Thus, sample size reduced to 32 in this latter model, and study site was no longer included as a random term. Change in management was not included because there were no significant differences in land use between years. We considered adding management itself as an explanatory variable as used in the first path model instead, which led to one alternative model to that presented in the results, where management replaced change in MST as the explanatory variable. Models were statistically not distinguishable (∆AICc ≤ 1.5), which was probably driven by the moderate correlation of change in MST and management (corr = −0.51, p = 0.003) in our study region. We further checked for the possible impact of the overall temperature trend over the decade (mean annual temperature 2009–2019) on species richness change as an alternative explanatory temperature variable to change in MST. Models were also not statistically different (∆AICc ≤ 1.5), but the explanatory value of the model got worse (R2 from 0.27 to 0.21), so that we concluded that change in MST appears to be the more suitable variable to explain the change in butterfly species richness. Therefore, we here present the (best) model, including change in MST, but discuss its uncertainty.

To identify shifts of communities and species, we included butterfly species with a total count of at least 10 individuals in each year and host plant species that occurred in at least 10 subplots per year. The exact magnitude of overall butterfly shifts as well as the number of included species varied according to the selected minimum number of individuals that a species needed to be included in the analysis, ranging from 17 included species (minimum 30 individuals) to 43 included species (minimum five individuals). For host plants, it varied between 44 included species (minimum 30 occupied subplots) and 85 included species (minimum five occupied subplots). The chosen minimum numbers of individuals/subplots reflected the overall change without excluding too many species (Appendix S1: Figures S4 and S5). Hence, we included 37 out of 71 butterfly species and 70 out of 128 host plant species.

We calculated butterfly community shifts by comparing the “community elevation score” between years, following Chen et al. (2009). To calculate the score, we assigned an average occurrence elevation to each species, using the abundance-weighted mean elevation of each species in 2009. Subsequently, we averaged the elevation values of all recorded species, weighted by abundance, to obtain the final score per site for each year. We analyzed the differences in community elevation score between years using a GAMM with the mgcv package version 1.8–38 (Wood, 2017). We fitted it with a smoothing function of site elevation, year as a factor (2009, 2019), and the interaction of site elevation and year, including study site as a random intercept effect.

To evaluate mean species range shifts, we used a linear mixed effect model within the lme4 package version 1.1–27.1 (Bates et al., 2015), with the abundance-weighted mean elevation per butterfly species as response variable and year as explanatory variable, including species as a random term. The same method was applied to check mean host plant species range shifts, replacing the response variable with the mean elevation of occupied subplots per host plant species. The difference in shifts between alpine and nonalpine species was calculated using a linear model with the shifts per species explained by the mean species elevation in 2009 and the factor alpine/nonalpine. A linear model with the elevations of all recorded individuals (butterflies)/occupied subplots (host plants) as response variable and year (continuous), species (factor), and their interaction as explanatory variables was used to test for shifts of single species (Figure 5). Exact shifts, corresponding standard errors and significance of shifts were derived post hoc by calculating estimated marginal means with the emmeans package version 1.7.2 (Lenth, 2022).

RESULTS

Interannual changes in temperature, butterfly, and host plant richness

In 2019, MAT on our study sites was on average 0.94°C warmer than in 2009. When comparing MST, it was 0.22°C warmer in 2019 than in 2009. In both cases, the difference in temperature between the two years changed significantly along the elevational gradient with a more pronounced temperature increase at higher elevations (GAMs, p ≤ 0.001, Figure 1 [dark lines] and Table 1). Interestingly, MST at elevations up to 1000 MASL was even cooler in 2019 than in 2009. The temperature differences between 2009 and 2019 also reflected the mean gradual temperature increase, which was even higher with 1.45°C in MAT and 1.23°C in MST over this decade. Equally, the gradual temperature increase over time changed significantly along the elevational gradient, with a stronger increase at high elevations (GAMs, p ≤ 0.001, Figure 1 [light lines] and Table 1).

Butterfly abundance did not change significantly along the elevational gradient but differed significantly between years (GAMM, elevation: year p ≤ 0.001, Table 1), being significantly lower in 2019 (GAMM, intercept between years p ≤ 0.001). The distribution of species richness along the elevational gradient changed from a linear decline in 2009 to a midelevation peak at ~1200 MASL in 2019 (GAMM, elevation: year p ≤ 0.05, Figure 2a and Table 1). Only in the lowlands (<1000 MASL) did species richness decline significantly when the elevational gradient was divided into three sections (lowland [644–999 m], midelevations [1000–1499 m], highlands [1500–2034 m]) and yearly species numbers were tested against each other within sections (t-tests, lowlands: t = 2.8776, df = 17.855, p = 0.0101, midelevations: t = −0.7474, df = 21.501, p = 0.4629, highlands: t = −0.6734, df = 18.863, p = 0.5089).

Details are in the caption following the image
Species richness of (a) butterflies and (b) host plants along the elevational gradient in 2009 (blue) and 2019 (orange) with 95% confidence intervals.

Altogether, we recorded 128 potential host plant species for the complete recorded butterfly species pool, of which we found 119 in 2009 and 122 in 2019. Host plant richness peaked at low and midelevations and declined above ~1250 MASL (Figure 2b). There was no significant difference in host plant richness between 2009 and 2019 along the elevational gradient (GAMM, difference between years p > 0.1, Table 1).

Drivers of butterfly richness

Path analysis indicated that MST, management, and host plant richness had a positive effect on butterfly richness (Figure 3a,b), with MST having the greatest overall effect. MST and host plant richness affected butterfly species richness directly, whereas extensive grassland management showed an indirect positive effect on butterfly richness by increasing host plant richness. Direct, negative effects of management on butterfly richness and positive effects of flower cover were not significant but included in the best-fitting path model (Fisher's C = 1.354, df = 2, p = 0.508, N = 64) that explained 63% of the variance in butterfly species richness. There was one alternative path model (∆AICc ≤ 1.5) where the nonsignificant relationship between MST and host plant diversity was discarded (dashed line). In this model, only the standard estimate of the relationship from management to host plant diversity changed from 0.47** (best model) to 0.62***.

Details are in the caption following the image
Drivers of butterfly richness along spatial and temporal gradients in National Park Berchtesgaden. (a) A priori hypothesized causal structure of the relationships among mean summer temperature (MST), management, larval (host plant diversity), and adult (flower cover) food resources and species richness of butterflies and (b) best fitted path model (based on the Akaike information criterion with a correction for small sample sizes, AICc, N = 64). The only alternative model (ΔAICc ≤ 1.5) had one path removed (dashed line), indicating an unstable path. (c) A similar causal structure was hypothesized for a path model that included the absolute differences in explaining variables between 2009 and 2019 (e.g., ΔMST = MST[2019]–MST[2009]), whereby change in management was excluded as there were no changes between years. (d) Again, the best-fitted path model is presented (N = 32). However, note that, although explanatory value was lower, the presented model was statistically not distinguishable (ΔAICc ≤ 1.5) from two other models replacing change in temperature with either overall temperature trend (mean 2009–2019) or management. The standardized path coefficients, their statistical significance (p ≤ 0.1 [.], p ≤ 0.05 [*], p ≤ 0.01 [**], p ≤ 0.001 [***]), and the conditional coefficients of determination (R2) are given. Paths with a significance of p ≤ 0.05 are presented in black and higher values in gray. Flower cover was (b) log-transformed in the first model and (d) square-root-transformed in the second one, prior to analysis.

To explicitly check the mechanisms that explain the difference (∆) in species richness between 2009 and 2019, we ran a second path model (Figure 3c,d). We only included the differences between years in all explaining and response variables per site. The best-fitting path model (Fisher's C = 0.194, df = 2, p = 0.907, N = 32) including these delta values explained 27% of the change in butterfly species richness between 2009 and 2019. In this model, the change in MST was the only variable significantly affecting the change in butterfly species richness between years. There was one alternative model (∆AICc ≤ 1.5) where the nonsignificant relationship between ∆MST and ∆ flower cover was discarded (single arrow model: ∆MST → ∆ butterfly species richness, standard estimate 0.43, p ≤ 0.05). However, it explained a lower proportion of the variance of the ∆ in species richness (R2 = 0.19). The same applied when replacing ∆MST with management or overall temperature trend over the decade (mean annual temperature 2009–2019): The path models were not statistically different from the one with ∆MST (∆AICc ≤ 1.5), but they explained less of the change in species richness (R2 = 0.23 or R2 = 0.21).

Upslope shifts of communities and species

The butterfly community elevation score, averaging the mean occurrence elevation of all detected butterfly species per site, showed a general upwards shift of species communities, as the score in 2019 was significantly lower than in 2009 (GAMM, elevation: year p ≤ 0.001, Figure 4 and Table 1). Averaging the difference in community elevation scores between years of each plot, communities shifted on average 51 ± 14 m (mean ± SD) upwards.

Details are in the caption following the image
Community elevation score of butterflies along the elevational gradient in 2009 (blue) and 2019 (orange) with 95% confidence intervals. A lower score in 2019 indicates that the community at a certain elevation is now hosting a higher share of species that were found at lower elevations in 2009.

On average, butterfly species moved 79 ± 19 m upwards based on a comparison of the difference in mean elevation between the two years (linear mixed effect model, df = 36, p ≤ 0.001). Shifting distances decreased significantly with increasing altitude (p ≤ 0.05), but in contrast, mountainous species (lower elevational limit ≥500 MASL according to Paolucci [2013]) shifted 139 ± 57 m more than other species (linear model, df = 34, p ≤ 0.05, Appendix S1: Figure S6). At the individual species level, 14 out of 37 selected butterfly species shifted significantly along the elevational gradient between 2009 and 2019 (estimated marginal means, post hoc to linear model [df = 6312, p ≤ 0.001]), 11 upwards and three downwards (Figure 5 left). Although most mountainous species shifted (eight out of 11), only about a quarter (six out of 26) of the species occurring in lowlands changed significantly. Host plants shifted with an average of 7 ± 8 m considerably less (linear mixed effect model, df = 69, p > 0.1) and also the proportion of shifting species was lower (Figure 5 right): Only six out of 70 host plants shifted significantly upwards or downwards (estimated marginal means, post hoc to linear model [df = 8079, p ≤ 0.001]), and most of them were grasses and belonged to either the Poaceae or Cyperaceae family.

Details are in the caption following the image
Abundance-weighted mean elevational shifts of butterfly (left) and host plant species (right) between 2009 and 2019 (vertical lines) with standard error (horizontal lines). The numbers in brackets give the number of individuals per species sampled in 2009/2019. Alpine butterfly species (lower elevational limit ≥500 m above sea level according to Paolucci (2013)) are shown in light and dark blue. Dark blue and black indicate significant shifts with respective significance levels shown next to the species names (left in each panel): * = p ≤ 0.05, ** = p ≤ 0.01, *** = p ≤ 0.001 (Wilcoxon rank-sum tests).

Considering the average increase of ~1.45°C in MAT over the last decade (Appendix S1: Figure S2) and a mean decrease in temperature of 0.52°C per 100 m of elevation (Appendix S1: Figure S3), the observed upslope shift of butterflies lagged well behind the expected shift of species required to track their temperature niche along the elevational gradient (~281 m). However, the reported gradual temperature increase over time also needs to be handled with care because the slope of the underlying linear model strongly depends on the years included in the model (e.g., MAT before 2009); here, we only considered a short time frame (2009–2019). Consideration of the absolute difference in MAT between 2009 and 2019 revealed that the expected shift was lower, about 179 m upwards, in accordance with an increase of ~0.93°C, which was closer to the observed shift rates.

DISCUSSION

Within a single decade, butterfly communities and species have changed along elevational gradients in ways consistent with predicted responses to temperature increases, whereas associated host plant communities did not change in the considered time frame.

Shifts in butterfly but not host plant richness

Although butterfly abundances were significantly lower in 2019 across the elevational gradient, the richness and composition of the regional species pool was very similar in both years. This contrasts with results from Spain, where significant overall declines in butterfly richness along elevational gradients over three decades were reported (Wilson et al., 2007) and signal the efficiency of the national park protection status in Berchtesgaden, in accordance with results from protected areas in the Italian Alps, where species richness even increased (Cerrato et al., 2019). Interestingly, we detected a clear upward shift of the maximum in butterfly species richness along the elevational gradient over time. Although in 2009 butterfly richness was highest in the lowlands, it peaked at midelevations in 2019. In our study, the hump-shaped pattern was mainly caused by the severe decline in species richness in the lowlands, whereby the majority of species were not lost but shifted to higher elevations, as generally expected under the gradually increasing temperatures detected in the region over the decade (Cerrato et al., 2019; Vitasse et al., 2021). Interestingly, within this period, the butterfly species decline in the lowlands was not compensated by species immigrating from even lower elevations or further south, similar to observations from Spain (Wilson et al., 2007). Such lowland biotic attrition has been described for the tropics, where no warmer habitats exist (Colwell et al., 2008). In temperate systems, it might indicate an unsuitable habitat matrix surrounding the national park, preventing butterfly species from following their climatic niches (Lenoir et al., 2020). Further, an intensification of anthropogenic influence might contribute to the lower species richness at lower elevations, as observed in several mountain ranges around the world (Forister et al., 2010; Gallou et al., 2017; Molina-Martínez et al., 2016). Our lowest study sites are located outside or at the outer margins of the national park. Thus, they might be more affected by land-use intensification (Viterbi et al., 2013) and, for example, higher amounts of insecticide use in the surrounding area (Brühl et al., 2021), even if local land use on the study sites stayed consistent.

In contrast to the shift in butterfly species richness patterns, we detected stable host plant richness along the elevational gradient. The hump-backed relationship between host plant species richness and elevation was also observed elsewhere (Gutiérrez et al., 2016) and matched the typical richness pattern observed for plants along elevational gradients in the Alps (Dainese & Poldini, 2012; Fontana et al., 2020; Hoiss et al., 2012). That there was no significant upward shift in richness might be due to the relatively time-limited study period. Plants react more slowly to temperature changes than mobile terrestrial insects (Vitasse et al., 2021). Over the entire last century, significant upward shifts in plant richness are visible (Lenoir et al., 2008), especially for mountainous and grassy species, which represent the majority of host plant species in our study. Also, an increase in plant richness at high elevations was observed, accelerating over time (Steinbauer et al., 2018). Thus, potential shifts might now become visible much faster than previously, which is supported by the slight tendency toward an upward shift of host plant richness already visible in our study over two years.

Temperature as main driver of butterfly richness

Butterfly species richness was directly driven by temperature (MST = mean summer temperature) and diversity of larval food resources (host plant richness), whereas host plant richness was driven by management, thereby indirectly also affecting butterfly richness. Adult food resource availability (flower cover) was not a significant predictor, but it did improve the path model, indicating its importance for butterfly species richness.

The strong positive correlation of butterfly richness with temperature is well known (Habel et al., 2021; Pellissier et al., 2013; Viterbi et al., 2013) as is the positive correlation of butterfly richness with (host) plant richness (Gutiérrez et al., 2016; Koch et al., 2013; Pellissier et al., 2013). However, in most cases, the two effects are examined separately without setting them in a common context. To analyze the drivers of butterfly richness in Britain, they were combined, finding a bigger impact of temperature in the form of growing degree-days than of host plant richness (Menéndez et al., 2007), corroborating our results. Nevertheless, land management was also not included in this study, which we found to be the main direct driver of host plant richness and, thus, indirectly driving butterfly richness. That a direct positive effect of extensive management on butterfly richness was found previously (Bubová et al., 2015; Jerrentrup et al., 2014) might be because (host) plant richness was not included as a possible intermediate factor. In addition, other factors not included in our analysis probably also exert an influence on butterfly richness, causing the remaining unexplained variance in our data. For example, structural diversity can positively affect butterfly richness (Viterbi et al., 2013), which we did not consider because data were not available for the first study year. Further, average annual precipitation has a positive effect on butterfly richness (Habel et al., 2021) but was not assessed at our study sites and is hard to predict at the site level for our region due to the complex topography. The latter might also be responsible for some temperature-correlated variance left unexplained, since microclimatic influences like topographic shading were not included, when modeling temperature for all sites, eventually causing some uncertainty in predictions.

According to the ascertained positive correlation between MST and butterfly richness, we would consequently expect an increase in butterfly richness with higher temperatures over the decade, and indeed, we identified the change in MST between years as the driver of change in butterfly richness between years as having the best explanatory power. Host plant richness did not explain the change; it also did not change between years. This finding is in accordance with other results explaining butterfly range shifts with temperature, but not with host plants (Wilson et al., 2005). However, in our study, the positive correlation between change in MST and change in butterfly richness resulted mainly from lower temperatures and similarly lower species numbers in the lowlands, whereas higher temperatures at higher elevations did not lead to an expected equally steep increase in species richness. The reason for this might be a time lag in the reaction of species richness to higher temperatures. This reaction might take more time than a decade; true abandonment of lower elevations or colonization of higher elevations would be required to observe changes in richness. This lag might also indicate that our results mirror delayed reactions to changes that had already occurred before the considered decade, suggesting a worse correlation between the observed changes in drivers and the reaction of species (richness). Further, we could not clearly separate the management and temperature effect, so that the negative response of species richness in the lowlands might also originate from management, possibly in interaction with temperature, as already shown in tropical mountains (Peters et al., 2019). Concerning the change in host plant richness, the missing link to change in MST in our path model seems to be a further indicator that the considered period might just have been too short to detect plant shifts. Thus, we speculate that butterfly richness changed to a certain degree with changing temperatures, whereas host plant richness had a bigger time lag in reacting because plants face a variety of challenges when trying to move upward, resulting in multiple types of lags (Alexander et al., 2018). Although less than in the lowlands (Bertrand et al., 2011), those lags could create an elevational gap between the peaks of butterfly and host plant richness in the future. Synthesizing other results from the Alps supports this assumption, as butterflies were the only order of studied organisms reacting in time to follow the upward-moving climatic niches (Vitasse et al., 2021).

Butterfly community and species upslope shifts

The community elevation score showed a clear upward shift of butterfly communities, as already reported from other regions in the world (Chen et al., 2011; Nieto-Sánchez et al., 2015; Wilson et al., 2007). The change was predominantly driven by losses of high-elevation species in the lowlands and range expansions of lowland species into midelevations, reflecting the observed changes in butterfly richness with losses in low elevations and an increase at midelevations. The slight inversion to higher community elevation scores in 2019 at the upper end of the elevational gradient might be a simple artifact because some alpine species monitored at the upper end of the elevational gradient might now have optimum ranges above the monitored elevational range. Alternatively, or additionally, temperature increases might also have increased the reproduction success and, thus, the abundances of alpine species. Because abundances are generally low at high altitudes, such an increase in high-altitude species abundances might outweigh the effect of upward-shifting lowland species.

Like the community elevation score, individual butterfly species moved considerably upward along the elevational gradient. The decreasing range shifts with increasing altitude are in line with the findings of previous studies (Mamantov et al., 2021; Roth et al., 2014), suggesting that species at higher altitudes might have certain traits slowing their shifts down (Cerrato et al., 2019; Mamantov et al., 2021; Rödder et al., 2021) or, conversely, that species at lower elevations are more pressured to react because they are closer to their upper thermal limits (Freeman et al., 2021). Interestingly, nevertheless, filtering for species restricted to higher elevations (only occurring at ≥500 MASL) shows that the mountainous species pool is reacting especially strongly (Figure 5 and Appendix S1: Figure S6), as also observed elsewhere (Cerrato et al., 2019; Molina-Martínez et al., 2016). It must be noted that the calculated significance and shifting distances of individual species were mainly abundance-driven, meaning that plots at lower elevations were not completely abandoned, but considerably more individuals were found in higher elevations, which is the reaction that could be expected over such a relatively short period. Some species did not show a significant response, possibly because of specific species traits. For example, widespread generalist species do not have a very specific distribution area that could shift visibly up- or downward.

Butterfly species in our study shifted stronger per decade than in other studies over 20–35 years (Forister et al., 2010; Molina-Martínez et al., 2016; Wilson et al., 2005). This might reflect the accelerating speed of climate change within the last decade. However, the shifts we found were still considerably below the predicted shifts to follow the climatic niches along the elevational gradient according to the observed temperature increase. Although we might have overestimated the predicted shifts owing to the limited considered time frame, a lag between observed and predicted shifts would also persist when assuming lower temperature increases in the Alps, as suggested by observations of the last 50 years (Vitasse et al., 2021). Lags were also reported in other studies (Chen et al., 2011; Devictor et al., 2012; Lenoir et al., 2020; Nieto-Sánchez et al., 2015) and appear more pronounced in temperate species communities, where species might partly benefit from higher temperatures (Freeman & Class Freeman, 2014). In temperate regions, species are still further from their upper thermal limits than in the tropics (Ghalambor, 2006), increasing the chances that species might just expand their ranges into higher elevations instead of moving their elevational range upward. Previous observations from the Alps strengthen this hypothesis, detecting stronger shifts in upper range limits than in the mean elevation of occurrence (Vitasse et al., 2021).

Nevertheless, the observed butterfly uphill shifts in mean elevation were still considerable (79 ± 19 m) and slightly above other recent observations from the Alps detecting ~40–50 m uphill shift/decade (Rödder et al., 2021; Roth et al., 2014; Vitasse et al., 2021). Thus, even if there might be an especially pronounced temporary response to the extraordinary hot year 2019 reflected in our data, a considerable species shift was found consistently throughout the Alps. At the same time, host plant species showed nearly no upward movement (7 ± 8 m), coinciding with the observations of the richness pattern and with other recent alpine studies showing similar range shifts for plants (Roth et al., 2014; Vitasse et al., 2021) when averaging shifts of alpine and herbaceous plants in the second study. The vast majority of host plants showed more uphill movement than grasses, which can be disregarded when looking at butterflies, because butterfly larvae feeding on grasses are mostly specialized at the family level, not on single plant species. Thus, moving away from a Poaceae or Cyperaceae species has no effect on them because they will find other related species. This fits together with our findings of a missing impact of host plant richness on changes in butterfly richness between years. Thus, when butterflies continue to follow their temperature niches uphill, this might form a gap between butterfly and host plant distributions, forcing them to adapt. A possible way for butterflies to avoid shifting away from their host plants is by adapting their phenology instead (Vitasse et al., 2021). The connection between host plant availability and this kind of adaptation has been demonstrated, since especially butterflies more dependent on certain larval host plants were observed to exhibit this strategy (Diamond et al., 2011). Otherwise, butterflies might also compensate for faster upward shifts broadening their larval diet breadth by including new plant species in their diet, as plants at higher elevations are less resistant to herbivores (Pellissier et al., 2012).

CONCLUSIONS

Within a single decade, we detected a clear upward shift of butterflies, concerning with maximum species richness, community composition, and the abundance-weighted mean elevation of individual species. While butterfly richness was explained by temperature and host plant richness, only a change in temperature correlated with changes in richness, though it was just driven by equally lower temperatures and richness in the lowlands, which might also be caused by management or land-use change in the surroundings. Since richness at higher elevations did not increase accordingly with increasing temperatures, we assume that there was a lag in reaction to higher temperatures, coinciding with our observations from individual species and community shifts, which also lagged behind the actual temperature change, despite showing considerable upward movements already. Because our study tracked community changes between two sampling years and not over an entire decade, we could not conclusively assess whether the observed butterfly community responses were reversible in cool years or ongoing responses to climate change. However, the general findings coincide well with the results from long-term studies over several decades, and warmer years are predicted to become more frequent in the near future. Thus, the reported responses in butterfly communities are likely to become the predominant ones. Importantly, larval host plants showed barely any shifts in the considered decade. We conclude that the broad distribution of host plants along the elevational gradient currently still allows fast butterfly responses to temperature changes. However, under ongoing climate change, a gap might open up between host plant and herbivore distributions, which could decelerate further butterfly shifts if they are not able to adapt to new host plants. Further, alpine species are facing continuous habitat loss, even without land-use change, as moving upward in a pyramid-shaped mountain range like the Alps goes hand in hand with a reduction in available habitat area. This underlines the importance of preserving traditional extensive pasture management practices in alpine regions since we showed that those support higher host plant and, thus, butterfly richness.

AUTHOR CONTRIBUTIONS

Alice Classen and Janika M. Kerner conceived the idea and developed the concept of the study. Janika M. Kerner conducted the butterfly monitoring in 2019, compiled the data, and conducted statistical analyses in an exchange of ideas with Fabienne Maihoff, Jochen Krauss, and Alice Classen. Fabienne Maihoff and Janika M. Kerner collected flower cover data and organized field work in 2019. Lukas Bofinger contributed to the vegetation assessment in 2019. Janika M. Kerner wrote the first draft of the manuscript, and all authors contributed to the final version.

ACKNOWLEDGMENTS

We sincerely thank the National Park Berchtesgaden for their support by granting access to the study area and providing housing. Special thanks to Ole Behling for coordinating logistics inside the park and to Annette Lotz and Daniela Kilian for cooperation in setting up the study. Further, we appreciated the cooperation with local farmers, who supported data collection on their pastures. We also thank Nikki Sauer for help in mounting voucher specimens of butterflies and Mirko Wölfling for his support in butterfly species identification. Equally, our gratitude goes to Roland Kaiser and Thomas Eberl from ENNACON, who conducted the vegetation survey. The study was funded by the Bavarian Ministry for Science and the Arts in the framework of the Bavarian Climate Research Network (bayklif). Open Access funding enabled and organized by Projekt DEAL.

    CONFLICT OF INTEREST

    The authors declare no conflict of interest.